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Abstract 

Three schemes, whose expressions are not too complex, are selected for 
the numerical integration of a system of stochastic differential equations 
in the Stratonovich interpretation: the integration methods of Heun, Mil- 
stein, and derivative- free Milstein. The strong (path-wise) convergence is 
studied for each method by comparing the final points after integrating 
with 2" and 2"~^ time steps. We also compare the time that the computer 
takes to carry out the integration with each scheme. Putting both things 
together, we conclude that, at least for our system, the Ifeun method is 
by far the best performing one. 



1 Introduction 

1.1 Stochastic differential equations 

Let us suppose that we are dealing with a system with deterministic dynamics, 
on which external noise is acting. That means that, in case we could remove the 
noise (by switching it off if possible, by isolating the system, etc) , its equations 
of motion would be purely deterministic: 

dXit)=fiX{t),t)dt. (1) 

Now, let us consider that Gaussian white noise is acting on the system. To 
include its effects, we add to (HJ a term g{X{t),t) dW{t): 

dX{t) = f{X{t), t) dt + g{X{t),t) dW{t), (2) 

where the functions of g are determined by some physical considerations, and 
the components of W are independent Wiener processes. A Wiener process P] is 
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a Gaussian stochastic process, almost surely continuous, where non-overlapping 
increments are independent, and with 



(i) 1^(0) = 0, 

(ii) {W{t)) = 0, 

(iii) VaT{W{t) ~ W{s)) = t-s. (3) 

1.2 Ito and Stratonovich interpretations 

Equation ([2|) can be written as 

X{t) = Xito) + f f{X{s),s) ds + f g{X{s),s) dT¥(s), (4) 

Jto Jto 

An issue arises when defining the integral over the Wiener processes. If we take 
a partition {to, . . . , tAr} of the interval [to, t] (where tjsf = t), and being ti a value 
between ti-i and ti, ti ^ [1 — A)ti_i + Xti, we would define the integral of any 
function $ as [5] 



/ <^{X{s),s)<lW{s) = lim V <^{X{n),n)[W{t,)-W{t,^^)], (5) 

where A = maxjt^ — The issue is that the limit above is different for 

different choices of A. The most common interpretations are the Ito scheme, in 
which the function is evaluated at the left-hand endpoint of each sub-interval 
(A = 0), and the Stratonovich scheme, in which the function is evaluated at 
the midpoint of each sub- interval (A = 1/2). We usually write the Stratonovich 
calculus with a circle before dW . 

The solutions to the stochastic differential equations are different for both 
interpretations, given that the stochastic integrals have a different value for each 
one. For example, the integral bellow in the Ito calculus is [T], 

f W{s) AW{s) = \[W{tf - W{tof] \{t to), (6) 
while for Stratonovich calculus, 

f W{s) o <lW{s) = \[W(tf - W{toY]. (7) 

Chosen one interpretation, one can always do the calculus in the other in- 
terpretation by modifying / with the Ito- Stratonovich drift correction formula 
[2 13] : the two equations bellow 



dX{t) = f{X{t),t)dt + g{X{t),t)dW{t) (8) 
dX{t) = l{X{t),t)dt + g{X{t),t)odW{t) (9) 

have the same solution if their drifts fulfil the relationship 
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l^{X{t),t)^MX{t),t)--gkj{X{t),t)—g,,{X{t),t). (10) 

When the equations (0) model a physical system, the "white noise" dW is 
actuaUy an ideahzation ol noise with a small correlation time: the autocorrela- 
tion function of dW{t) is not really a delta function, but a sharply peaked one; 
i.e., its correlation time is small, but positive. Therefore, the stochastic function 
dW{t) is actually not singular, and the stochastic differential equations ([2]) have 
a well-defined solution. In the limit of correlation time of dW{t) going to zero, 
such solution tends to the integration of the equations with the Stratonovich 
interpretation (see more details and a discussion in chapter IX. 5 of [4], Discus- 
sion of the Ito- Stratonovich dilemma). That is why we will use the Stratonovich 
interpretation. 

2 Integration methods 

We will consider three light-weight (in terms of the complication of the formulae) 
numerical integration methods leading to the Stratonovich interpretation. They 
all require using a sample of the discretized Wiener process at each integration 
step. Let h be the length of the time step used for the numerical integration. 
Therefore, according to the properties of the Wiener process ([3]), we need a 
sample of independent random variables ^i{tn) at each of the values of the time 
{tn} in the discretization, distributed as VhAf{0, 1) = A/'(0, h), where Af{a, a"^) 
stands for the normal distribution with mean a and variance ct^. 

2.1 Heun method 

One of the simplest discretization schemes leading to the Stratonovich interpre- 
tation is the Heun method [51 IHl [H [Z]- This is a predictor-corrector method: 
given the value of at a time tn of the discretization, we first obtain the 
predictors, or supporting values, with the Euler integration scheme 

Xi{tn+l) = Xi{tn) + fi{x{tn),tn) k + gtj {x{tn) , tn) {tn) , (H) 

where tn+i — tn + h. Then, we obtain X(tn+i) as 

Xi{tn+l) = Xi{tn) + ^[fi{x{tn),tn) + Mx(tn+l),tn+l)] h 

+ ^ [gi3{x{tn),tn) + gij{x{tn+l),tn+l)] (,j{tn)- (12) 

Note that we are using Einstein's tensor convention of summing over repeated 
indices. 

2.2 Milstein scheme 

The Milstein scheme requires the use of the derivatives of the diffusion coeffi- 
cients. In the general case, it reads [T] 

Xi{tn+l) ^ X^{tn) + fih + gij£,j + gij^^S^ J{j,k), (13) 
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where all the terms at the right hand side of ([T^ are taken at time t„, and J{j.k) 
is a multiple Stratonovich integral 

J[m = dW^,(s2)) dW^fc(si) (14) 

(when using the Ito interpretation, the corresponding multiple Ito integral is 
defined in the same way as in but with the Ito interpretation, and written 
as%fc).)^ 

These integrals cannot be easily expressed in terms of the increments and 
of the components of the Wiener process. Nevertheless, we can consider the 
particular case of diagonal noise, where each component of the system has its 
own, independent, noise, i.e.: there are as many independent Wiener processes 
(number of components of W) as the number of variables in the system (the 
number of components of X), each component Xi of X is affected only by the 
corresponding component Wi of the Wiener process, and the diagonal diffusion 
coefficient ga depends only on Xi and maybe on time (but not on the other 
components oi X). In such a case, gij in (fT3)) is not equal to zero only where 
I = j, gik only where i = k, and dga/dxi only where i = I. Putting all things 
together, the last term in (IT^ is equal to 

dgu J 

(it is clear that we are not summing over i in the expression above , given that i 
is an external index in (jl3p ; we will omit this remark wherever it would be easy 
to figure out over which repeated indices we are not summing up.) Regarding 
the last factor in the product above, (IT4l) becomes for j = fc = i. 



J(M) = J m{si)-W,{t^)]odW,{si) 

-t 



W,{si)om,{si)-W,{tn) odW,(si) 



= \[W^{t^+lf - W,{t^f] - M^.(t„)[M^,(i„+i) - W,{t^)] 

= \mitn+i)-W,{tn)]''. (15) 

where ([7]) has been used at the third equality. Therefore, for the numerical 
integration, we set [I] 
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2 

As a result, for diagonal noise the Milstein method dT5|) becomes, in the Stratonovich 
interpretation. 



Xi{tn+l) = Xi{tn) + ft h + gu 6 + ^ ^"^^ ^^"^"^ 



It is also possible to simplify the general expression ([13)) for the case of 
commutative noise jlj, a case slightly more general than our case of diagonal 
noise. 



4 



2.3 Derivative-free Milstein method 



The Milstein method above requires the analytic specification of the first deriva- 
tive of the diffusion term. This can be sometimes inefficient, either because the 
analytic expression of the derivative is highly complex, or because we are try- 
ing many different functions for g, etc. In such cases, we can use a numerical 
approximation for the derivative of g for use in the Milstein scheme. 

Here, we will only give the expression for diagonal noise. The formulae for 
the Ito interpretation can be seen on equations (1.3) and (1.4) of chapter 11 in 
[1]. First of all, the supporting values are obtained as 

Xi ^ Xi + fi h + guVh- (18) 
Then, starting from (1171) for the Stratonovich interpretation, we get 

Xi{tn+i) = Xi{tn) -f fih + gue_., + --^[gu{x,tn) - gu] (19) 

2v " 



3 Check of strong convergence 

For each method, we are interested in checking its degree of strong convergence, 
or convergence over the path. That means that, for any single realisation W^(i) 
of the Wiener process on a given time interval [t^, T\ - that must be chosen as the 
biggest time interval that we will use among all our simulations -, and starting 
with given initial conditions X[t^) for the system, the final point X{T) obtained 
with a good integration scheme, using a number of time steps big enough, must 
be close to the real solution of the differential equation, at time T (for that 
particular realisation of the Wiener process VF(t) and for this particular initial 
condition X(to)). Of course, when selecting any other realisation of the Wiener 
process and any other initial condition, the final point at time T obtained by 
the integration method must be close to the one of the real solution as well. 

Of course, we cannot know the value of X{T) obtained with the real solution 
but, still, we can test whether the integration method gives something close 
enough to it, by studying the self-consistence. If we have chosen a small number 
of integration steps (a big value for the time step /i), the obtained X{T) will 
not be reliable and, when repeating the integration with more time steps (but 
the same Wiener process and the same initial conditions) , we will obtain a new 
X{T) which will considerably differ from the previous one. On the other hand, 
we can be reasonably sure that we have already chosen a number of integration 
steps big enough, and therefore we have obtained a reliable value of X{T) if, 
after repeating the integration once more with a considerably bigger number 
of time steps (e.g., twice the former number of time steps, with same Wiener 
process and same initial conditions), the new obtained X{T) stays close to the 
value previously obtained. 

We are using here this method of self-consistency via the Brownian tree: 
we will start with the maximum number of timesteps 2^, for a chosen natural 
number iV, and we will then progressively reduce the number of integration 
steps by a half each time, i.e., we will integrate with 2^, 2^"-'^, 2^~^, etc, time 
steps, until the desired minimum power of 2 (that can be, if desired, as little as 
2°, i.e., a single time step). Given the discretised Wiener process {'Ci(2", j), j = 
0, . . . , 2" — 1} for 2" time steps , we obtain the discretisation of the same Wiener 
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process for half the number of time steps by summing up the two members of 
each couple, i.e. 

6(2"-\j)=e.(2", 2j)+6(2",2j- + l), j = 0,...,2"-i-l. (20) 

Sometimes, one may first take a signal of Gaussian white noise ^ with mean 
zero and variance 1 (i.e., ( is distributed as A/'(0, 1)), so the integration routine 
gets the increments ^ of the Wiener process from the 1-correlated signal C as 
^ = \/hC^, where h is the time step that is used at that moment to integrate 
the differential equation (see introduction to section [5]) . Then, when halving 
the total number of time steps (and therefore doubling the time step to 2/i), the 
new signal C still has to be 1-correlated, so we have to make 



C.(2"-i,j-) = ^[C(2", 2j)+C.(2", 2j + l)], j- = 0,...,2"-i-l, (21) 

so that (^i(2"^^) is l~correlated as well as Ci(2"). Also, that way, given that 
C(2") = v^C(2") and C(2""^) = V2hC{T^~^), equation ^ yields 

For a desired precision for the i-th component Xi in the stochastic differ- 
ential equation, we will consider that 2"~^ is a number of steps big enough if 
\xi{T, 2") — Xi{T, 2"~^)| < A;. In order to check that such accuracy holds for 
different realisations of the noise, we will repeat the comparison of the 2"- and 
2"~ ^-integrations with different brownian paths {Wk |fc = 1, . . . ,K}: we call 
Xi{t,2",k) to the trajectory with the Wiener process Wk- We will always use 
the same initial conditions {xi(to)}- 

As a result, we have for each component i and each exponent n — 1 a set of K 
values consisting of the differences between the state at the final time T with 2" 
integration steps minus the state with 2"'"^ integration steps: each value corre- 
sponds to each realisation of the Wiener process. To be clear: we have for each 
component i and each exponent n a set Diff(i, n) = {I?fe(*, n) |fc = 1, . . . , K}, 
where Vkiiju) = Xi(T, 2"+^,fc) — Xi{T,2^ , k). From each set I?fc(i,?T,), we can 
assess the reliability of the integration with 2" time steps. Normally, if n is not 
too small, the mean value of those differences (2?) will be around 0, otherwise 
we can say that there is a preferred sign in the values of the differences V^ii, n) 
and, therefore, a noticeable systematic error at the numerical integration. Then, 
the quality of the numerical integration is better the smaller the absolute values 
of the differences \'Dk{i,n)\, i.e., the narrower the distribution, and this being 
peaked around 0. A way of measuring this is to obtain the mean and some 
central momenta, 

^l^{^,n)^{ {Vk{t,n) - (2?fe(^,r^)) )™ ) , (22) 

of the distribution. Given that, for a number of integration steps big enough, 
the absolute values of the differences \T>k{i, n)| should be as small as possible, we 
can say that the reliability is better the smaller (in absolute value) the mean and 
the central momenta. Of course, a number of time steps that is acceptable for 
a given integration method (Heun, Milstein, etc) will be in general insufficient 
or, contrarily, higher than necessary, for another integration method. 
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4 The system 



4.1 General description of the model 

We are using here a dynamical system as a model for paleoclimate. Our system 
consists of three variables: the volume of the ice V, the CO2 concentration C, 
and the variable D that is related to the ocean's temperature and circulation. 
The oscillatory astronomical forcing acting on the system can be approximated 
by two functions [8]: Il{t) carrying the effect of the precession of the Earth's axis, 
and E{t) carrying the effect of the obliquity. Both functions are approximated 
by a sum of harmonic oscillations, ^ sin(a;it + 9i), where the parameters a^, 
Wi and 9i are arranged as rows in two data files (a file for precession and another 
one for obliquity) , and we sum for each function the number of oscillations that 
we like according to the desired precision and computation time. 

Nine parameters appear in the equations that govern the deterministic sys- 
tem: an offset Vr in the coupling of the variable D to the ice volume, a coupling 
Ve of the ice volume to the obliquity, a coupling Vp of the ice volume to the 
precession, an offset Vq in the coupling to the ice volume, a coupling Cp of the 
CO2 to the precession, the relaxation times of each variable ry, tc, td, and a 
parameter ap playing a role in the coupling to the variable V . 

We will add to each equation a Gaussian white noise variable: r]v{t), r]c{t), 
r]D{t) to the variables V, C, D respectively. The functions rj are independent, 
with zero mean and variance one (the actual variances of the noises will be 
included in the functions that multiply the ?7's): 



{V^{t)) = 

{v^{t)v,{t')) = 5,,5{t~t'), (23) 

(these functions 77 correspond to the Wiener processes AW in ([5]).) After in- 
serting the noise, we will also have to take into account the variances of the 
noises, and maybe some more extra parameters in the functions that couple the 
variables V , C, D to the noises (i.e., the functions that multiply the noises). 
Putting all things together, the equations of the system read: 



dV 
'dt 


- — [^y(V)- 
Tv 


^R] + 


dC 
"d^ 


- — {C + V- 


-Id- 
2 


dD 
'df 


td 


-{V- 



where the function R is defined as: 



(24) 



R = afl(e^° - 1) + (1 - afl,)i?o 

Ro = VpR + VeE + ^C+^D + Vo, (25) 



and the potentials are: 
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0.04 



V 

1 3 

-X^ - X. 



(26) 



As stated before, the variance of each noise is taken into the function that 
multiplies it. In the simplest case of additive noise, we just make gi = Ui, where 
af is the variance of the noise acting on the component i. 

4.2 Manual corrections 

It is worth mentioning that we are performing two corrections "by hand" at each 
integration step. On the first hand, in order to prevent 1^ < 0, which would be 
devoid of physical meaning, we shall make V = 0.001 by hand whenever we end 
up with V < 0.001 after an integration step. On the other hand, in an attempt 
to reduce computing divergences, we will bound fv (V) a.nd make it actually be 
^v{V) = -min{4, 0.04/^}. 

4.3 Optimization 

When coding, we have of course tried to avoid redundant calculations as much 
as possible. This specially includes the figures that just depend on time (and 
not on the state variables V, C, D). Given that, in most cases, we will be 
propagating many particles, we will first obtain once and for all the values that 
depend only on time and are therefore the same for all the particles. We start by 
generating the array {t} of all the times at which the integration steps will take 
place, and then we calculate at all such times: the precession 11, the obliquity E, 
and also the combination (Vp H + Ve E + Vq) appearing in (^5)1 . This will imply 
passing these three extra parameters to all the functions (instead of passing just 
the parameter time) and adds complexity to the code, but it saves invaluable 
time when running the programme, specially when the number of particles is 
high. 

4.4 Particular choices in this paper 

For the tests that will follow in the next sections, we took 50 terms for the 
precession and 20 terms for the obliquity. Given that we are comparing three 
methods of numerical integration, we do not want to stay in the simple case of 
additive noise, but we want to test the different methods in the more general 
case of multiplicative noise. We therefore use a simple kind of multiplicative 
noise: 



The values of the parameters used here for the comparison of the numerical 
integration schemes are the following (one unit of time corresponding to 1000 
years): 



9v{V,t) 
9c{C,t) 
9D{D,t) 



<TvV 

<JcC 
<JdD. 



(27) 
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Vt = 0.9, 
Cp = 
TV = 19, 
ap = 0.3 

2 2 



Tc = 10, Tb = 1 




Ve = 0.14, = 0.21, Vo = 0.82 



= 0.001. 



(28) 



The initial conditions used are ^(0) = 0.33, C(0) = 0.5, £'(0) = 0. 

5 Result of the comparison of the integration 
schemes 

5.1 Computing scenario 

The code was written in C and compiled with the Intel compiler (ice) in order 
to take advantage of vectorisation. The compilation commands were 

ice && ice file.c -L/opt/gsl/lib -Igsl -Igslcblas 
-I/opt/gsl/include -o file 

(we use the first command ice to switch to 64 bits.) 

The CPU is Intel (R) Xeon(R) X5450 @ 3 . OOGHz running on SUSE Linux 
11.0 x86_64. 

5.2 Generation of the same trajectories from the three 
integration schemes 

Trajectories for three particles were generated using the three integration schemes, 
with the same Wiener process for the three schemes. We integrated up to time 
equal to 2000, using 10^ time steps per trajectory. 

The time used to generate the data file for each integration scheme on the 
computing scenario described above was: 

• Hcun scheme: 0.335 seconds. 

• Milstein scheme: 0.304 seconds. 

• Derivative-free Milstein scheme: 0.357 seconds. 

Nevertheless, such small times are not reliable for comparing the different in- 
tegration schemes in terms of time expenses, given that such execution times 
vary from one execution to another of the same executable file, depending on 
the other jobs in the whole computer. We will compare the time expenses in 
Section 15.31 

In order to check the agreement between the different integration schemes, 
we have plotted, in figures [1] and [21 a two-by-two superposition of the same 
trajectories generated with the different integration schemes. For clarity, we 
have only plotted the second half, from time 1000 to time 2000. The space 
between consecutive points has been intentionally made relatively large, in order 
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Figure 1: Superposition of the trajectories generated by the Heun and Milstein 
integration schemes. The same number in the key corresponds to the same 
trajectory. 
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1- heun 
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1\ n -ft \ 
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Figure 2: Superposition of tire trajectories generated by the Heun and 
derivative-free Milstein integration schemes. The same number in the key cor- 
responds to the same trajectory. 
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to be able to see the overlapping in the same trajectory obtained from two 
different integration schemes. 

We can see in figures [T] and [2] that all the trajectories overlap perfectly 
with the corresponding same trajectory obtained with the other integration 
scheme. This says that the three integration methods are reliable for the number 
of integration steps stated above (10^). Now, we will make in section [^751 a 
quantitative study of accuracy versus time expense of each integration scheme. 

5.3 Rate of strong convergence for the three integration 
methods 

To compare the rate of strong convergence of the three integration schemes, 
we followed the method described in section [31 We integrated 10'* particles 
from time to time 400 with each integration scheme. The maximum number 
of time steps was 2^^, and the minimum 2^ (the biggest power of 2, as we 
will see, that is small enough to cause divergences in the numerical integration 
for the three schemes). The same Wiener processes were used for the three 
integration schemes. For the sample of differences, at the final time 400, between 
the integrations with 2" and 2"~^ time steps, we obtain the mean and central 
momenta (l22t . for 2 < m < 4. 

The time used to generate the data file for each integration scheme on the 
computing scenario described above was: 

• Heun scheme: real 4m44.192s, user 4m29.721s, sys 0m4.396s. 

• Milstein scheme: real 3m45.887s, user 3m41.430s, sys 0m4.008s. 

• Derivative-free Milstein: real 4m21.324s, user 4m9.640s, sys 0m4.296s. 

Given that the most time consuming task of the programme is the numerical 
integrations, we can say from the figures above that, for the same number of 
time steps, the Milstein scheme is the fastest method; the next one would be 
the derivative-free Milstein taking around 1.13 times longer than the former, 
and the slowest one is the Heun scheme taking around 1.22 times longer than 
the Milstein scheme (beware, though, that these ratios may change if we use 
another computing environment). Given that these ratios are close to 1 (i.e., 
the time expense of the integration schemes is roughly the same for the three), 
we shall choose the most accurate of our integration schemes. 

In order to compare the accuracies, we have rearranged the tables to display 
together the same relevant variables derived after integrating with the three dif- 
ferent methods. The means and central momenta explained above are displayed 
in Tables [THU The columns are named after the system's variable (V, C or D), 
and after the integration method: suffix "-h" for Heun, "-m" for Milstein, and 
"-df" for derivative-free Milstein. The names of the rows correspond to "the 
small power of 2", i.e., n means that we are considering the distribution of the 
differences, at the final time 400, between the integrations with 2"+-'^ and 2" 
time steps: small absolute values in the mean and central momenta indicate 
that the integration with 2" time steps is in principle reliable. Also, for the sake 
of visual clarity, all the means and central momenta have been multiplied by 
10^ before displaying them in the tables. 
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;) 


V-h 


V-m 


V-df 


C-h 


C-in 


C-df 


D-h 


D-ni 


D-df 


15 


47 


485 


-1485 


-41 


-336 


1133 


20 


464 


-1183 


14 


-31 


1560 


-1301 


13 


-1142 


1060 


-4 


1585 


-638 


13 


136 


2097 


-1490 


-156 


-1870 


1034 


-76 


-187 


-967 


12 


226 


4931 


-862 


-134 


-4070 


970 


117 


2832 


797 


11 


1511 


10451 


2723 


-1136 


-8673 


-2193 


1222 


5354 


1869 


10 


2.4 e5 


nan 


nan 


-2.3 e5 


nan 


nan 


-1.2 e5 


nan 


nan 


9 


nan 


nan 


nan 


nan 


nan 


nan 


nan 


nan 


nan 



Table 1: Mean values times 10^; for the different variables, integration methods, 
and time steps. See text for an explanation on the notation. 



n 


V-h 


V-m 


V-df 


C-h 


C-m 


C-df 


D-h 


D-m 


D-df 


15 


41 


487 


710 


24 


441 


710 


5 


475 


1101 


14 


31 


596 


393 


22 


464 


421 


4 


1872 


1573 


13 


56 


923 


837 


36 


851 


595 


321 


3693 


2499 


12 


484 


1661 


1090 


255 


1515 


1099 


1886 


3307 


4130 


11 


1008 


3181 


2806 


990 


2796 


2269 


2997 


9920 


7689 


10 


1.3 o5 


nan 


nan 


9.6 c4 


nan 


nan 


3.1 e5 


nan 


nan 


9 


nan 


nan 


nan 


nan 


nan 


nan 


nan 


nan 


nan 



Table 2: Second central momenta times 10^; for the different variables, integra- 
tion methods, and time steps. 



n 


V-h 


V-m 


V-df 


C-h 


C-m 


C-df 


D-h 


D-m 


D-df 


15 


4 


-364 


-80 


-3 


322 


-198 





524 


-1503 


14 


-10 


314 


-63 


4 


-17 


57 





3120 


-691 


13 


9 


-121 


-272 


-6 


122 


125 


-542 


-3701 


-1759 


12 


-56 


29 


-407 


54 


161 


534 


-47 


3415 


2488 


11 


93 


72 


-81 


43 


288 


224 


1347 


5900 


3729 


10 


-33385 


nan 


nan 


17172 


nan 


nan 


-36116 


nan 


nan 


9 


nan 


nan 


nan 


nan 


nan 


nan 


nan 


nan 


nan 



Table 3: Third central momenta times 10^; for the different variables, integration 
methods, and time steps. 



n 


V-h 


V-m 


V-df 


C-h 


C-m 


C-df 


D-h 


D-m 


D-df 


15 


5 


544 


457 


2 


408 


632 





815 


2253 


14 


5 


344 


98 


2 


205 


132 





5454 


3934 


13 


7 


371 


247 


3 


372 


130 


956 


9851 


6762 


12 


352 


844 


374 


61 


758 


522 


5938 


6916 


9426 


11 


360 


1353 


2232 


386 


1077 


1206 


7854 


23677 


19313 


10 


1.1 e5 


nan 


nan 


6.6 e4 


nan 


nan 


8.4 e5 


nan 


nan 


9 


nan 


nan 


nan 


nan 


nan 


nan 


nan 


nan 


nan 



Table 4: Fourth central momenta times 10*^; for the different variables, integra- 
tion methods, and time steps. 
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We can see that the numerical integration yields divergencies when done 
with a small number of time steps: in our case, they appear when we descend 
to 2^^ integration steps for the two Milstein methods, and when we descend to 
2^ steps for the Heun method. As we have stated in Section H?^ we have tried 
to remove the divergences coming from the equations, so we should think that 
the divergences found are due in principle only to the integration using too few, 
too big, time steps. 

5.4 Choice of the numerical integration method 

If we check a given row (corresponding a given number of integration steps) in 
Tables [T}]31 we see that the mean and central momenta for the Heun scheme 
are always smaller (in absolute value) compared to the figures for the other two 
schemes (with a couple of exceptions in Table [3]). Not only smaller, but the fig- 
ures for the Heun scheme are almost always (except for some rows corresponding 
to a small number of time steps) at least one order of magnitude smaller than 
the figures for the two Milstein schemes. If we were too strict, one might object 
that the Heun scheme is (just slightly) more expensive in terms of computing 
time: however, in most cases, the figure for the Heun scheme in a given row 
of the tables is smaller than the figure for the other schemes displayed one row 
above (which corresponds to an integration with double number of time steps) , 
so the slightly bigger time expense of running the Heun scheme is largely made 
up for by the accuracy of the method. 

As a result, we can state that the Heun method is the best performing as 
we get much better accuracy, for a given computation time, compared to the 
Milstein methods. The better performance of the Heun method was to some 
extend expectable, given that, in the case of additive noise, the Milstein scheme 
reduces to the basic Euler scheme (as the derivative in (jl7p vanishes), whereas 
Heun's scheme is always (also for additive noise, and even for no noise at all) a 
second-order predictor-corrector method. 

Also, the idea of Heun's method is easy to understand, and the method is 
easy to code (in particular, it does not use any derivative). All this makes the 
Heun's scheme a very suitable and attractive method in our opinion. 
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